La accidentalidad vial es un fenómeno que se está convirtiendo en un problema de salud pública dado que toca directamente al ser humano y detrás de cada cifra hay un drama familiar, por esta razón las entidades gubernamentales tienen que ser más ambiciosos en materia de campañas, entendiendo que hay una gran inversión desde el orden nacional para resolver esta problemática, pero no es suficiente, por lo cual se debe insistir más en el tema de seguridad vial desde lo local. Al respecto, la alcaldía de Medellín a venido recolectando datos sobre accidentalidad vehicular en el municipio desde el año 2014 con el fin de describir el comportamiento y los lugares más afectados por este problema. En este trabajo se utiliza información de accidentes de 2014 a 2018 para describir estadisticamente los accidentes presentados y tratar de predecir accidentes presentados en las diferentes epocas del año y sitios del municipio.
La Estadística se ocupa del manejo de la información que pueda ser cuantificada. Implica esto la descripción de conjuntos de datos y la inferencia a partir de la información recolectada de un fenómeno de interés. La función principal de la estadística abarca: Resumir, Simplificar, Comparar, Relacionar, Proyectar. Constantemente se generan datos que se alojan en los sistemas de información y bases de datos públicas y privadas. La Alcaldía de Medellín cuenta con una página web donde alojan bases de datos públicas en diferentes temas sobre la ciudad. El objetivo de este trabajo es utilizar la Estadística y el desarrollo computacional para analizar los datos de incidentes georreferenciados de 2014 a 2018 de la ciudad de Medellín para predecir accidentes del año 2019 con el fin de descubrir factores que se relacionan con los accidentes y ayudar de alguna manera a las entidades públicas a tomar decisiones en temas de seguridad vial.
En el enlace se encuentra alojado el aplicativo web VehApp que permite visualizar la accidentalidad y pronosticar Enlace: Aplicacion web VehApp
En el enlace se encuentra alojado el Github principal que contiene todos los archivos y scripts utilizados para la construcción e implementación de este trabajo Enlace: Github principal
En inteligencia artificial un problema se define como la existencia de discrepancias entre un estado actual y un estado objetivo. La resolución recurrente y especializada de ciertos tipos de problemas permite identificarlos en forma de tareas. Para abordar este problema se usa la metodología para la Analítica de datos (ASUM-DM) ya que establece una forma probada para llevar un problema a una solución de analítica. ASUM-DM de IBM se basa en la metodología CRISP-DM y hace enfasis en las nuevas prácticas en la ciencia de datos como el uso de volumenes de datos muy grandes, la incorporación de análisis de texto, en el modelado predictivo y la automatización de algunos procesos. El ciclo de elaboración se basa en 9 principales pasos:
Más de un millón de personas muere cada año en accidentes automovilísticos, Estos incidentes ocurren por una cantidad de razones, pero el error humano está en el centro. Los conductores manejan más rápido de lo que deberían en ciertas condiciones atmosféricas, toman decisiones equivocadas y no logran anticipar riesgos potenciales. El conocimiento de las variables involucradas en los accidentes de tránsito pueden facilitar la formulación de políticas específicas para reducir la creciente accidentalidad vial. Asi, el objetivo de este trabajo es describir las características de tiempo, lugar, vehiculo y demás variables que puedan inferir en los accidentes de tránsito de la ciudad de Medellín entre 2014 y 2018, buscando aportar elementos al conocimiento de este problema que sean útiles para quienes allí formulan las políticas de prevención y control de morbilidad y mortalidad.
La Alcaldía de Medellín en el portal GeoMedellín suministra las bases de datos anuales para incidentes de tránsito registrados por la Secretaria de Movilidad. Los datos se encuentran en este enlace.. Para el entendimiento del conjunto de datos se define el accidente de transito como: “evento, generalmente involuntario, generado al menos por un vehiculo en movimiento, que causa daños a personas y bienes involucrados en él e igualmente afecta la normal circulación de los vehiculos que se movilizan por la via o vias comprendidas en el lugar o dentro de la zona de influencia del hecho”
La estructura de las variables que contiene los datos es la siguiente:
| Variable | Descripción | Tipo | Detalle |
|---|---|---|---|
| X | Coordenada X (longitud) de la ubicación del accidente | float | - |
| Y | Coordenada Y (latitud) de la ubicación del accidente) | float | - |
| OBJECTID | Identificación del registro | integer | - |
| RADICADO | Identificación única del accidente ante la Secretaría de Movilidad | string | - |
| FECHA | Fecha de la ocurrencia del accidente | string | - |
| HORA | Hora aproximada del accidente | datetime | - |
| DIA | Número día dentro del mes del accidente | integer | - |
| PERIODO | Año de la ocurrencia del accidente | integer | - |
| CLASE | Tipo de accidente | String | Clases: atropello, caída del ocupante, choque, incendio, volcamiento, otro |
| DIRECCION | Dirección detallada de la ubicación del accidente | String | - |
| DIRECCION_ENC | Dirección encasillada de la ubicación del accidente | string | Formato único de direcciones de la Alcaldía de Medellín |
| CBML | Identificación única del lote más cercano a la ubicación del accidente | string | Acrónimo de comuna, barrio, manzana, lote |
| TIPO_GEOCOD | Tipo de ubicación catastral | string | Formato de la oficina de catastro de la Alcaldía de Medellín |
| GRAVEDAD | Gravedad del accidente | string | Clases: herido, muerto, solo daños |
| BARRIO | Barrio donde ocurrió el accidente | string | 340 niveles dentro de la variable |
| COMUNA | Comuna donde ocurrió el accidente | string | 85 niveles dentro de la variable |
| DISENO | Tipo de diseño en donde ocurrió el accidente | string | Clases: ciclo ruta, glorieta, intersección, lote o predio, paso a nivel, paso elevado, paso inferior, pontón, puente, tramo de vía, túnel, vía peatonal |
| DIA_NOMBRE | Nombre del día de la semana del accidente | string | - |
| MES | Número de mes dentro del año del accidente | integer | - |
| MES_NOMBRE | Nombre del mes del accidente | string | - |
| X_MAGNAMED | - | - | - |
| Y_MAGNAMED | - | - | - |
| LONGITUD | Coordenada longitud de la ubicación del accidente | float | - |
| LATITUD | Coordenada latitud de la ubicación del accidente | float | - |
El conjunto de datos llamado medellin cuenta con 228.693 registros y 24 variables que corresponden a accidentes desde el 2014 a 2018. Para el inicio de la preparación de ellos se analizará que variables no hay necesidad de tener en cuenta, cuantas variables hay que transformar y que otras se pueden agregar. Para esto se hace uso de Tidyverse en el lenguaje R para la limpieza de datos y otros paquetes necesarios.
Desde el portal de GeoMedellín se descargan datos de accidentalidad por año desde 2014 a 2018 en archivos csv, se ingresan a una carpeta del computador y se hace lectura de cada uno de ellos luego se juntan para crear solo un conjunto de datos que abarca los 5 años.
# Cargar libreria
library(data.table);library(purrr);library(dplyr)
library(plotly);library(tidyr);library(stringr)
library(lubridate);library(sf);library(knitr);library(kableExtra)
## lectura de datos
# lista archivos
setwd("Datos")
archivos_datos <- list.files(pattern = ".csv", include.dirs = T, recursive = F)
# lectura de los archivos (2014 a 2018)
datos <- map(archivos_datos, fread, sep = ",", encoding = "UTF-8", colClasses = "c")
# juntar por filas para crear un solo conjunto de datos
medellin <- bind_rows(datos)
# ver datos
kable(head(medellin), align = "ccccc") %>%
kable_styling() %>%
scroll_box(width = "100%", height = "300px")| X | Y | OBJECTID | RADICADO | FECHA | HORA | DIA | PERIODO | CLASE | DIRECCION | DIRECCION_ENC | CBML | TIPO_GEOCOD | GRAVEDAD | BARRIO | COMUNA | DISENO | DIA_NOMBRE | MES | MES_NOMBRE | X_MAGNAMED | Y_MAGNAMED | LONGITUD | LATITUD |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| -75.6027257778941 | 6.21901553077807 | 211279 | 1423828 | 2014/01/01 00:00:00+00 | 02:20:00 | 01 | 2014 | Choque | CR 80 CL 8 | CR 080 008 000 00000 | 1611 | Malla vial aproximada: CR 81-7 | HERIDO | Loma de los Bernal | Belén | Tramo de via | MIÉRCOLES | 01 | 831190.08 | 1179689.5 | -75.60272578 | 6.21901553 | |
| -75.5681847049393 | 6.26000899508826 | 211280 | 1423839 | 2014/01/01 00:00:00+00 | 12:50:00 | 01 | 2014 | Atropello | CR 53 CL 61 | CR 053 061 000 00000 | 1003 | Malla vial | HERIDO | Jesús Nazareno | La Candelaria | Tramo de via | MIÉRCOLES | 01 | 835013.63 | 1184223.79 | -75.56818471 | 6.260009 | |
| -75.5499431117657 | 6.26476471445961 | 211281 | 1423840 | 2014/01/01 00:00:00+00 | 01:00:00 | 01 | 2014 | Atropello | CR 39 CL 70 | CR 039 070 000 00000 | 0308 | Malla vial | HERIDO | Manrique Oriental | Manrique | Tramo de via | MIÉRCOLES | 01 | 837032.71 | 1184749.86 | -75.54994311 | 6.26476471 | |
| -75.6076107982279 | 6.23432736551888 | 211282 | 1423849 | 2014/01/01 00:00:00+00 | 12:37:00 | 01 | 2014 | Atropello | CL 32 CR 84 | CL 032 084 000 00000 | 1617 | Malla vial | HERIDO | Las Mercedes | Belén | Tramo de via | MIÉRCOLES | 01 | 830649.47 | 1181383.24 | -75.6076108 | 6.23432737 | |
| -75.5796854467144 | 6.29996777075499 | 211283 | 1423890 | 2014/01/01 00:00:00+00 | 10:40:00 | 01 | 2014 | Caída de Ocupante | CR 80 CL 101 | CR 080 101 000 00000 | 0603 | Malla vial | HERIDO | Doce de Octubre No.2 | Doce de Octubre | Tramo de via | MIÉRCOLES | 01 | 833740.78 | 1188643.8 | -75.57968545 | 6.29996777 | |
| -75.5554293007481 | 6.28449892412236 | 211284 | 1423892 | 2014/01/01 00:00:00+00 | 04:00:00 | 01 | 2014 | Choque | CR 48 CL 93 | CR 048 093 000 00000 | 0401 | Malla vial | HERIDO | Berlin | Aranjuez | Interseccion | MIÉRCOLES | 01 | 836425.43 | 1186932.72 | -75.5554293 | 6.28449892 |
Con la descripción de las variables que contienen los datos se identifica las siguientes correcciones que se deben hacer a los datos:
Convertir la FECHA en formato YYYY-MM-DD
Eliminar datos vacios en la variable CLASE, quitar tildes y unificar nombres de niveles
Los datos vacios en la variable DISENO se cambian por OTRO, quitar tildes y unificar nombres de niveles
Quitar tildes a nombres de dias
# Organizar fecha
medellin$FECHA <- as.Date(medellin$FECHA) # Formato YYYY-MM-DD
# Eliminar datos vacios de la variable CLASE y correccion de niveles
medellin <- medellin[-which(medellin$CLASE == ""),]
# Correccion tildes
medellin$CLASE <- iconv(medellin$CLASE, from = "UTF-8", to = "ASCII//TRANSLIT")
# Quitar "DE" (nivel "Caída de ocupante" variable CLASE)
medellin$CLASE <- gsub(pattern = "Caida de Ocupante", replacement = "Caida Ocupante", x = medellin$CLASE)
#Cambiar nivel vacio por otro en variable DISENO
medellin$DISENO <- ifelse(medellin$DISENO == "","Otro",medellin$DISENO)
# Quitar tildes en niveles de la variable DISENO
medellin$DISENO <- iconv(medellin$DISENO, from = "UTF-8",to="ASCII//TRANSLIT") #Quitar tildes
# Quitar tildes a nombres de dias en variable DIA
medellin$DIA_NOMBRE <- iconv(medellin$DIA_NOMBRE, from = "UTF-8", to = "ASCII//TRANSLIT")
# ver estructura variables
str(medellin)## Classes 'data.table' and 'data.frame': 228687 obs. of 24 variables:
## $ X : chr "-75.6027257778941" "-75.5681847049393" "-75.5499431117657" "-75.6076107982279" ...
## $ Y : chr "6.21901553077807" "6.26000899508826" "6.26476471445961" "6.23432736551888" ...
## $ OBJECTID : chr "211279" "211280" "211281" "211282" ...
## $ RADICADO : chr "1423828" "1423839" "1423840" "1423849" ...
## $ FECHA : Date, format: "2014-01-01" "2014-01-01" ...
## $ HORA : chr "02:20:00" "12:50:00" "01:00:00" "12:37:00" ...
## $ DIA : chr "01" "01" "01" "01" ...
## $ PERIODO : chr "2014" "2014" "2014" "2014" ...
## $ CLASE : chr "Choque" "Atropello" "Atropello" "Atropello" ...
## $ DIRECCION : chr "CR 80 CL 8" "CR 53 CL 61" "CR 39 CL 70" "CL 32 CR 84" ...
## $ DIRECCION_ENC: chr "CR 080 008 000 00000" "CR 053 061 000 00000" "CR 039 070 000 00000" "CL 032 084 000 00000" ...
## $ CBML : chr "1611" "1003" "0308" "1617" ...
## $ TIPO_GEOCOD : chr "Malla vial aproximada: CR 81-7" "Malla vial" "Malla vial" "Malla vial" ...
## $ GRAVEDAD : chr "HERIDO" "HERIDO" "HERIDO" "HERIDO" ...
## $ BARRIO : chr "Loma de los Bernal" "Jesús Nazareno" "Manrique Oriental" "Las Mercedes" ...
## $ COMUNA : chr "Belén" "La Candelaria" "Manrique" "Belén" ...
## $ DISENO : chr "Tramo de via" "Tramo de via" "Tramo de via" "Tramo de via" ...
## $ DIA_NOMBRE : chr "MIERCOLES" "MIERCOLES" "MIERCOLES" "MIERCOLES" ...
## $ MES : chr "01" "01" "01" "01" ...
## $ MES_NOMBRE : chr "" "" "" "" ...
## $ X_MAGNAMED : chr "831190.08" "835013.63" "837032.71" "830649.47" ...
## $ Y_MAGNAMED : chr "1179689.5" "1184223.79" "1184749.86" "1181383.24" ...
## $ LONGITUD : chr "-75.60272578" "-75.56818471" "-75.54994311" "-75.6076108" ...
## $ LATITUD : chr "6.21901553" "6.260009" "6.26476471" "6.23432737" ...
## - attr(*, ".internal.selfref")=<externalptr>
La variable barrio y comuna es muy importante para el objetivo, por eso se hace un análisis completo de como están estructurados sus niveles y las reformas que se le deben hacer. Se encontró que la variable COMUNA tiene nombres erroneos y además contiene 19.736 registros vacios. En la variable BARRIO se tienen 19.766 registros vacios y algunos de sus niveles contienen numeros que no indican bien el nombre de la ubicación. Para la solución de este problema se usó la variable CBML que es una identificación única del lote más cercano al accidente conformado por un Acrónimo de comuna, barrio, manzana, lote, con esta variable se puede extraer un registro con los 4 primeros números para obtener la comuna y el barrio. La oficina de catastro cuenta con un shapefile donde se registran verdaderamente los nombres de barrios y comunas identificados con su codigo que en los datos de accidentes serían esos 4 primeros digitos del CBML, así que se hace un cruce entre los datos de accidentes y los datos de catastro uniendo por el código. Como observación la variable CBML cuenta con 19.757 datos vacios y 185 registros en los cuales no cuenta con un registro CBML claro, que ayude a identificar el barrio y la comuna.
LATITUD Y LONGITUD dado que X y Y representan lo mismo, RADICADO, CBML, TIPO GEOCOD, X_MAGNAMED, Y_MAGNAMED, DIRECCION, DIRECCION_ENC y MES_NOMBRE también se eliminan y se realiza un ajuste en el nombre de las variables. El nuevo conjunto de datos queda compuesto de 208.745 registros y 15 variables
#==========================
# ARREGLOS COMUNA Y BARRIO
#===========================
setwd("Datos")
catastro <- read_sf("Limite_Barrio_Vereda_Catastral.shp")
# Crear un codigo desde CBML que deje solo comuna y barrio
medellin <- mutate(medellin, CB = str_sub(CBML,1,4))
# ingresar las varaibles de identificacion de lugar dadas por catastro en el archivo de datos
medellin <- inner_join(medellin, select(catastro, CODIGO, NOMBRE_COM, NOMBRE_BAR),
by = c("CB" = "CODIGO"))
# Quitar duplicados por el inner_join. Pasamos de 228687 registros a 208745 registros
medellin <- medellin[!duplicated(medellin[,c("OBJECTID")])]
# Quitar tildes y volver a mayuscula los nombres de barrios
medellin$NOMBRE_BAR <- iconv(toupper(medellin$NOMBRE_BAR), from = "UTF-8", to = "ASCII//TRANSLIT")
# Elegimos las columnas que necesarias
medellin <- medellin[,c("X","Y","OBJECTID","FECHA","HORA","DIA","PERIODO","CLASE","GRAVEDAD",
"DISENO","DIA_NOMBRE","MES","CB","NOMBRE_COM","NOMBRE_BAR")]
# Cambiemos algunos nombres de las variables
names(medellin)[c(1,2,14,15)] <- c("LONGITUD","LATITUD","COMUNA","BARRIO")
# Asi quedan los datos
head(medellin)En las nuevas variables a ingresar se considera días especiales que corresponden a fechas que se caracterizan por ser grandes eventos en los que la comunidad participa de manera completa o mayoritaria y que de alguna manera afecta la movilidad, para esto se crea un archivo plano csv que tiene el registro de cuales han sido fechas especiales desde 2014 a 2018 y el tipo de festividad, para esta variable se tuvo en cuenta las siguientes fechas especiales.
Otra variable nueva a ingresar se crea a partir de la hora del accidente. En la ciudad de Medellín se considera como hora a pico al tiempo en que la mayoría de los habitantes se dirigen al trabajo, estudio o de regreso a casa, asi el intervalo de tiempo que se toma como hora pico es de 5am a 8am y de 5pm a 8pm. Cabe resaltar en la estructuración de la variable hora no está especificado a que parte del día ocurrió el accidente, es decir, con la hora no se puede saber si ocurrió en la mañana, en la tarde o en la noche.
Otras dos variables nuevas a ingresar corresponden al estado del clima:
HUMEDAD_REL: Humedad relativa - Es una relación, expresada en porcentaje, de la cantidad de humedad atmosférica presente relativa a la cantidad que estaría presente si el aire estuviera saturado, esta relación es determinante para que se produzca evaporación o precipitación (lluvia)PLUVIO: Pluviometria [mm] - Mide la cantidad de precipitaciones (lluvias) caídas en un lugar durante un tiempo determinadoSe debe advertir que la variable HUMEDAD_REL tiene registros desde 2014-12-04 a 2018-12-01 y la variable PLUVIO desde 2014-01-01 a 2018-12-01
Los nuevos datos quedan conformados por 208.745 observaciones y 19 variables. Para finalizar con el arreglo de los datos se realiza una estructuracion en el tipo de variable necesaria para los análisis.
#==============================
# Ingreso de fechas especiales
#==============================
setwd("Datos")
# Lectura registros de fechas especiales desde el 2014 a 2019
especiales <- read.csv("Otros/Fechas_especiales.csv", sep = ";", header = T)
# Convertir en formato YYYY-MM -DD
especiales$FECHA <- as.Date(especiales$FECHA, format = "%d/%m/%Y")
# unir las fechas especiales a medellin, los que nos son fechas especiales los pone como NA
medellin <- merge(x = medellin, y = especiales, by = "FECHA", all.x = T)
medellin$ESPECIAL <- ifelse(is.na(medellin$ESPECIAL),"NO","SI")
# Variable HORA_PICO
medellin$HORA <- hms(medellin$HORA)
medellin$HORA_PICO <- ifelse(hour(medellin$HORA) %in% c(5,6,7), "SI","NO")
# Variable HUMEDAD_REL y PLUVO
# Humedad Relativa
HUMEDAD <- read.csv('Otros/Datos IDEAM/HumedadRelativa.csv.csv')
HUMEDAD <- data.frame(FECHA=as.Date(HUMEDAD$Fecha),HUMEDAD_REL=HUMEDAD$Valor)
# Pluviometria
PLUV <- read.csv('Otros/Datos IDEAM/Pluviometria.csv.csv')
PLUV <- data.frame(FECHA=as.Date(PLUV$Fecha),PLUVIO = PLUV$Valor)
PLUV <- PLUV[!duplicated(PLUV[,c("FECHA")]),]
medellin <- merge(x = medellin, y = HUMEDAD, by = "FECHA", all.x = TRUE)
medellin <- merge(x = medellin, y = PLUV, by = "FECHA", all.x = TRUE)
# Tipo variables
medellin$LONGITUD <- as.numeric(medellin$LONGITUD)
medellin$LATITUD <- as.numeric(medellin$LATITUD)
medellin$DIA <- as.integer(medellin$DIA)
medellin$PERIODO <- as.integer(medellin$PERIODO)
medellin$CLASE <- as.factor(medellin$CLASE)
medellin$GRAVEDAD <- as.factor(medellin$GRAVEDAD)
medellin$DISENO <- as.factor(medellin$DISENO)
medellin$DIA_NOMBRE <- as.factor(medellin$DIA_NOMBRE)
medellin$MES <- as.integer(medellin$MES)
medellin$COMUNA <- as.factor(medellin$COMUNA)
medellin$BARRIO <- as.factor(medellin$BARRIO)
medellin$ESPECIAL <- as.factor(medellin$ESPECIAL)
medellin$HORA_PICO <- as.factor(medellin$HORA_PICO)
# Visualizar datos
head(medellin)Y su nueva estructura de variables es la siguiente:
| Variable | Descripción | Tipo | Detalle |
|---|---|---|---|
| FECHA | Fecha ocurrencia accidente | Date | - |
| LONGITUD | Coordenada longitud ubicación accidente | integer | - |
| LATITUD | Coordenada latitud ubicación accidente | integer | - |
| OBJECTID | Identificación única del registro | string | - |
| HORA | Hora aproximada del accidente | datetime | - |
| PERIODO | Año de la ocurrencia del accidente | integer | - |
| CLASE | Tipo de accidente | FACTOR | 7 niveles: Atropello, Caida Ocupante, Choque, Choque y Atropello, Incendio, Volcamiento, Otro |
| GRAVEDAD | Tipo de gravedad del accidente | FACTOR | 3 niveles: HERIDO, MUERTO, SOLO DAÑO |
| DISENO | Tipo de diseño en donde ocurrió el accidente | Factor | 13 niveles: Ciclor Ruta, Glorieta, Interseccion, Lote o Predio, Paso a nivel, Paso Elevado, Paso Inferior, Ponton, Puente, Tramo de via, Tunel, Via peatonal, Otro |
| DIA NOMBRE | Nombre del dia de la semana del accidente | Factor | 7 niveles: LUNES, MARTES, MIERCOLES, JUEVES, VIERNES, SABADO, DOMINGO |
| MES | Número de mes dentro del año del accidente | integer | - |
| CB | Identificación del lote más cercano a la ubicación del accidente | string | Acrónimo de Comuna y barrio |
| COMUNA | Comuna donde ocurrió el accidente | Factor | 21 niveles |
| BARRIO | Barrio donde ocurrió el accidente | Factor | 307 niveles |
| ESPECIAL | Fecha especial | Factor | niveles: SI, NO |
| HORA_PICO | Hora pico en la ciudad (5am - 8am, 5pm-8pm) | Factor | niveles: SI, NO |
| HUMEDAD_REL | Porcentaje Humedad relativa | Numeric | - |
| PLUVO | Cantidad de precipitaciones [mm] | Numeric | - |
En los últimos años ha habido un acelerado incremento del parque automotor circulante en las ciudades del país, lo que genera problemas de movilidad, tales como la congestión y aumento de la accidentalidad. De acuerdo con Censo de Población y Vivienda 2018 realizado por el DANE en el país hay 48.258.494 habitantes de los cuales se estima el 30.3% tenía hasta Junio 20 de 2019 algún tipo de vehiculo registrado en el parque automotor donde el 17.2% de esa población tenía principalmente motocicletas, el 7.1% automoviles y el 3.3% camionetas.
En medellín, según cálculos del Registro Único Nacional de Tránsito (Runt) el parque automotor del Valle de Aburrá viene en crecimiento continuo durante los últimos años. Solo entre 2016 y 2019 el número de vehiculos se incrementó en más de 210.000 unidades, arrojando un total de 1.674.554 vehiculos en la ciudad. Desde 2014 a 2018 en la ciudad se presentaron 228.693 accidentes, en los cuales ocasionaron la muerte a 1.374 personas, 126.222 resultaron heridas y además en 101.097 accidentes dejaron solo daños materiales. El siguiente análisis descriptivo se realiza teniendo en cuenta el objetivo de este trabajo y cabe resaltar que en la preparación de los datos se eliminaron 19.948 registros repartidos en los 5 años debido a duplicados o ausencia del dato.
require(ggplot2)
require(ggpubr)
#============= Análisis Temporal
#===============
# Dias
#===============
# Registros por día
medellin_dia <- medellin %>% group_by(DIA_NOMBRE) %>% summarize(total_accidentes = n())
# Ordenar de Lunes a Domingo
medellin_dia$DIA_NOMBRE <- ordered(medellin_dia$DIA_NOMBRE,
levels = c("LUNES","MARTES","MIERCOLES","JUEVES","VIERNES",
"SABADO","DOMINGO"))
grafico_1 <- ggplot(data = medellin_dia, aes(x = DIA_NOMBRE, y = total_accidentes)) +
geom_bar(stat = "identity", alpha = 0.65, fill = "grey20", color = "black") +
geom_text(aes(y = total_accidentes, label = total_accidentes),
position = position_dodge(width = 0.8), size = 2.8, vjust = 0.5, hjust = -0.1, col = "gray10") +
xlab("Días") +
ylab("Total Accidentes") +
ggtitle("Total Accidentes por día de la semana") +
theme_minimal() +
ylim(c(0,40000)) +
coord_flip()
#=================
# Mes
#=================
nombre_mes = c("Ene","Feb","Mar","Abr","May","Jun","Jul","Ago","Sept","Oct","Nov","Dic")
med_mes <- table(medellin$MES)
names(med_mes) <- nombre_mes
med_mes <- as.data.frame(med_mes)
grafico_2 <- ggplot(data = med_mes, aes(x = Var1, y = Freq)) +
geom_bar(stat = "identity", alpha = 0.7, fill = "grey20", col = "black") +
geom_text(aes(y = Freq, label = Freq),
position = position_dodge(width = 0.7), size = 2, vjust = -0.5, hjust = 1.5, col = "gray10") +
xlab("Mes") +
ylab("Total Accidentes") +
ggtitle("Total Accidentes por mes") +
ylim(c(0,20000)) +
theme_minimal()
#=================
# Anio
#=================
med_Anio <- table(medellin$PERIODO) %>% as.data.frame()
grafico_3 <- ggplot(data = med_Anio, aes(x = Var1, y = Freq)) +
geom_bar(stat = "identity", alpha = 0.65, fill = "grey20", col = "black") +
geom_text(aes(y = Freq, label = Freq),
position = position_dodge(width = 0.7), size = 3.5, vjust = 1.2, hjust = 0.5, col = "gray10") +
xlab("Año") +
ylab("Total Accidentes") +
ggtitle("Total Accidentes por año") +
theme_minimal()
#===================
# Semana
#===================
med_sem <- medellin %>% group_by(semana=week(medellin$FECHA)) %>%
summarise(Total_Accidentes = n())
grafico_4 <- ggplot(data = med_sem, aes(x = semana, y = Total_Accidentes)) +
geom_bar(stat = "identity", alpha = 0.6, fill = "grey20", col = "black") +
xlab("Año") +
ylab("Total Accidentes") +
ggtitle("Total Accidentes por Semana") +
theme_minimal()
# Imprimir los 4 graficos
ggarrange(grafico_3, grafico_2, grafico_4,grafico_1, ncol = 2, nrow = 2)#======================
# Serie temporal diaria
#=======================
med_diaria <- medellin %>% group_by(FECHA) %>% summarize(total_accidentes = n())
med_diaria$Anio <- year(med_diaria$FECHA)
med_diaria$dia <- day(med_diaria$FECHA)
med_diaria$mes <- month(med_diaria$FECHA)
plot_ly(data = med_diaria, x = ~FECHA, y = ~total_accidentes,
type = "scatter", mode = "lines", split = ~Anio,
line = list(width = 1)) %>%
layout(title = 'Registros de accidentalidad diarios 2014-2018',
xaxis = list(title = "Día"),
yaxis = list(title = "Número de accidentes"))El mes con menos promedio de accidentes diarios es Enero, con un promedio de 97 accidentes, mientras que Agosto presenta el mayor promedio de accidentes diarios con 121, esto teniendo en cuenta los 5 años de registros. El domingo es el día de la semana que menos accidentes se presentan con un promedio de 75 mientras que los demás días están alrededor de 100 y 140 accidentes en promedio. El día que presenta mayor variabilidad en cantidad de accidentes es el viernes, presentando varios datos atipicos por encima y por debajo.
# Gráfico de lineas Promedio diario mensual por anio
aggregate(total_accidentes~Anio*mes, data = med_diaria, FUN = mean) %>%
plot_ly(x = ~mes,
y = ~total_accidentes, type = "scatter", mode = "lines",
split = ~Anio, line = list(width = 1.5)) %>%
layout(title = 'Promedio diario mensual por año',
xaxis = list(title = "Mes"),
yaxis = list(title = "Número de accidentes"))# Boxplot por dia de la semana
med_diaria$semana <- weekdays(med_diaria$FECHA)
med_diaria$semana <- factor(med_diaria$semana,
levels = c("lunes","martes","miércoles","jueves","viernes","sábado","domingo"))
plot_ly(data = med_diaria, x = ~semana, y = ~total_accidentes,
type = "box", marker = list(color = 'grey24')) %>%
layout(title = 'Accidentes por día de la semana',
xaxis = list(title = "Día de la semana"),
yaxis = list(title = "Accidentes"))Se analiza la variable barrio y comuna, La Candelaria es la Comuna que mas presenta accidentes entre 2014 y 2018 seguida por Poblado. El Barrio donde más se presenta accidentes es justamente es el barrio La Candelaria, también se puede notar del mapa que la zona centro es la más afectada por los accidentes y las vias principales.
# Accidentes por comuna
medellin_comuna <- medellin %>% group_by(COMUNA) %>%summarize(accidentes = n())
ggplot(data = medellin_comuna, aes(x = reorder(COMUNA,+accidentes), y = accidentes)) +
geom_bar(stat = "identity", position = "dodge", fill = "grey24", color = "black", alpha = 0.6) +
geom_text(aes(y = accidentes, label = accidentes),
position = position_dodge(width = 0.7), size = 3.5, vjust = 0.5, hjust = -0.1, col = "black") +
xlab("Comuna") +
ylab("Total accidentes") +
ggtitle("Total accidentes por comuna 2014-2018") +
ylim(c(0,50000)) +
theme_minimal() +
coord_flip()#================================================
# MAPA
#================================================
require(leaflet)
# Mapa para todos los barrios
nueva_base <- medellin %>% filter(PERIODO >= 2014 & PERIODO <= 2018) %>%
group_by(CB) %>%
summarise(accidentes = n()) %>%
ungroup()
nuevo_mapa <- inner_join(catastro, nueva_base, by = c("CODIGO" = "CB"))
# Paleta de color
mypal <- colorNumeric(palette = c("#000000","#280100","#3D0201","#630201","#890100","#B00100","#DD0100","#F50201",
"#FF5F5E","#FF7A79","#FF9796","#FEB1B0","#FDC9C8", "#FFE5E4"), domain = nuevo_mapa$accidentes, reverse = T)
# Crear mapa
leaflet() %>% addPolygons(data = nuevo_mapa, color = "#0A0A0A", opacity = 0.6, weight = 1, fillColor = ~mypal(nuevo_mapa$accidentes),
fillOpacity = 0.6, label = ~NOMBRE_BAR,
highlightOptions = highlightOptions(color = "black", weight = 3, bringToFront = T, opacity = 1),
popup = paste("Barrio: ", nuevo_mapa$NOMBRE_BAR, "<br>", "Accidentes: ", nuevo_mapa$accidentes, "<br>")) %>%
addProviderTiles(providers$OpenStreetMap) %>%
addLegend(position = "bottomright", pal = mypal, values = nuevo_mapa$accidentes, title = "Accidentes", opacity = 0.6)Solo el 0.6% de los accidentes ha dejado muertos el otro 99.4% ocasionan Heridos o solo daños. El 88% de los accidentes ocurren en un Trammo de vía o intersección y la clase de accidente más frecuente es el Choque, además los choques son los que generan más daños y el atropello es el que genera más muertos junto a los choques, la frecuencia de heridos se nota en todos los tipos de accidentes.
# ------------ Gravedad del accidente
med_gravedad <- table(medellin$GRAVEDAD) %>% prop.table() %>% as.data.frame
med_gravedad$Var1 <- ordered(med_gravedad$Var1, levels = c("HERIDO","SOLO DAÑOS","MUERTO"))
grafico_5 <- ggplot(data = med_gravedad, aes(x = Var1, y = Freq)) +
geom_bar(stat = "identity", alpha = 0.7, fill = "firebrick", color = "black") +
geom_text(aes(y = Freq, label = paste(round(Freq*100,1),"%")),
position = position_dodge(width = 0.8), size = 4, vjust = -0.5, hjust = 0.5, col = "gray10") +
xlab("Gravedad") +
ylab("Frecuencia (%)") +
ggtitle("Total Accidentes por gravedad") +
theme_minimal()
#-------- Diseño de la via en el accidente
med_diseno <- table(medellin$DISENO) %>% prop.table() %>% as.data.frame
med_diseno$Var1 <- ordered(med_diseno$Var1, levels = c("Via peatonal","Tunel","Ponton","Paso a Nivel","Paso Inferior","Puente","Ciclo Ruta",
"Paso Elevado","Otro","Glorieta","Lote o Predio","Interseccion","Tramo de via"))
grafico_6 <- ggplot(data = med_diseno, aes(x = Var1, y = Freq)) +
geom_bar(stat = "identity", alpha = 0.7, fill = "firebrick", color = "black") +
geom_text(aes(y = Freq, label = paste(round(Freq*100,2),"%")),
position = position_dodge(width = 0.8), size = 3.2, vjust = 0.5, hjust = -0.2, col = "gray10") +
xlab("Diseño") +
ylab("Frecuencia (%)") +
ggtitle("Total Accidentes por diseño de la via") +
theme_minimal() +
coord_flip() +
ylim(0,1)
#------------ Clase Accidente
med_clase <- table(medellin$CLASE) %>% prop.table() %>% as.data.frame
med_clase$Var1 <- ordered(med_clase$Var1, levels = c("Choque","Otro","Atropello","Caida Ocupante","Volcamiento","Incendio","Choque y Atropello"))
grafico_7 <- ggplot(data = med_clase, aes(x = Var1, y = Freq)) +
geom_bar(stat = "identity", alpha = 0.7, fill = "firebrick", color = "black") +
geom_text(aes(y = Freq, label = paste(round(Freq*100,3),"%")),
position = position_dodge(width = 0.8), size = 4, vjust = -0.5, hjust = 0.5, col = "gray10") +
xlab("Clase") +
ylab("Frecuencia relativa (Accidentes)") +
ggtitle("Total Accidentes por Clase") +
theme_minimal() +
ylim(0,0.8)
# Imprimir los 3 graficos
ggarrange(grafico_5, grafico_6, ncol = 2, nrow = 1)# Clase vs gravedad
med_clase_gravedad <- table(medellin$CLASE, medellin$GRAVEDAD) %>% prop.table() %>% as.data.frame()
med_clase_gravedad$Var1 <- ordered(med_clase_gravedad$Var1, levels = c("Choque y Atropello","Incendio","Volcamiento","Caida Ocupante","Atropello","Otro","Choque"))
names(med_clase_gravedad) <- c("Var1","Gravedad","Freq")
ggplot(data = med_clase_gravedad, aes(x = Var1, y = Freq, fill = Gravedad)) +
geom_bar(stat = "identity", position = "dodge") +
geom_text(aes(y = Freq, label = paste(round(Freq*100,1),"%")),
position = position_dodge(width = 0.8), size = 3, vjust = 0.5, hjust = -0.2, col = "gray10") +
theme_minimal() +
ylim(0,0.5) +
xlab("Frecuencia (%)") +
ylab("Clase de Accidente") +
coord_flip()Con el fin de agrupar barrios según sus caracteristicas de accidentalidad se usa k-means para buscar los barrios mas cercanos a dichas caracteristicas. Para este trabajo las variables que caracterizan a los barrios son las siguientes:
Para obtener el número k optimo que reduce la distancia entre cluster, se realiza un gráfico de codos en el que se nota que el número k optimo es 4 u 5, asi que elegimos un \(k = 5\) dado la cantidad de barrios de la ciudad.
# Seleccionar las variables que se usan
agrupamiento <- medellin[,c("FECHA","BARRIO","COMUNA","PERIODO","CB","ESPECIAL","HORA_PICO",
"GRAVEDAD")]
# Agrupar las variables por barrijo
data_agrupamiento <- agrupamiento %>% filter(PERIODO >= 2014 & PERIODO <= 2018) %>%
group_by(CB) %>%
summarise(total_accidentes = n(),
Acc_Especiales = sum(ESPECIAL == "SI"),
Acc_Hpico = sum(HORA_PICO == "SI"),
Acc_Muerto = sum(GRAVEDAD == "MUERTO")) %>%
ungroup()
# Grafico de codos para k optimo
wssplot <- function(data, nc = 15, seed = 1234){
wss <- (nrow(data)-1)*sum(apply(data,2,var))
for (i in 2:nc){
set.seed(seed)
wss[i]<- sum(kmeans(data, centers = i)$withinss)}
plot(1:nc, wss, type = "b", xlab = "Number of Clusters",
ylab = "Within groups sum of squares", main = "Busqueda de numero de cluster optimo")}
wssplot(data_agrupamiento[,-1], nc = 10)Ahora, teniendo en cuenta que el numero de grupos o cluster optimo es 5, entonces se procede a usar kmeans de R para agrupar los barrios en 5 grupos y se realiza un grafico con sus dos primeras componentes principales para visualizar el agrupamiento. En la gráfica se puede notar que la dimension 1 que corresponde a el total de accidentes por barrio explica la variabilidad en un 79.5%.. Los datos de agrupamiento se guardan en un archivo csv llamado Agrupamiento_barrios.csv
require(factoextra)
require(tidyverse)
require(ggplot2)
# Ejecucion modelo para agrupamiento con kmenas
set.seed(1017239703)
barrios <- data_agrupamiento$CB
k.means.fit <- kmeans(data_agrupamiento[,-1], 5)
d <- k.means.fit$cluster
data_agrupamiento$Grupos <- d
# Visualizacion de grupos con sus dos primeras componentes principales
agrup_grafica <- data_agrupamiento[,-1]
row.names(agrup_grafica) <- barrios
set.seed(123456789)
fviz_cluster(object = k.means.fit, data = agrup_grafica,
ellipse.type = "convex", repel = F, show.clust.cent = FALSE,
labelsize = 8) +
labs(title = "Clasificación Jerárquica Barrios Medellín",
subtitle = "K Nearest Neighbor") +
theme_bw() +
theme(legend.position = "bottom")Ahora, para visualizar la agrupacion en el mapa haremos uso del shape file de catastro y de leaflet para dibujar el mapa y con colores representar los grupos seleccionados. En el mapa se puede notar claramente como los barrios del centro son los mas afectados por la accidentalidad y los barrios cercanos a las vias rapidas de la ciudad también son muy afectados. La idea del mapa es poder servir como gestor para implementación de estrategias que puedan reducir la accidentalidad de una forma más precisa, como ejemplo en la comuna La Candelaria que es la mas afectada por los accidentes y la mortalidad pero que cuenta con 20 barrios en los que se encuentra el Centro administrativo o el barrio Las Palmes que hacen parte del grupo de los mas bajos accidentes, pero se encuentra, Guayaquil, La Candelaria, Barrio Colón, Perpetuo Socorro y San Benito que hacen parte del selecto grupo de los de mayor accidentalidad en Medellín.
setwd("Datos")
data_agrupamiento <- read.csv(file = "Otros/Agrupamiento_barrios.csv", sep = ",", encoding = "UTF-8", colClasses = "character")
# Mapa para todos los barrios
dato_mapa <- inner_join(catastro, data_agrupamiento, by = c("CODIGO" = "CB"))
# Definicion de colores
thalia <- c("#FF0000", "#CCFF00", "#00FF66", "#0066FF", "#CC00FF")
dato_mapa$colores <- ifelse(dato_mapa$Grupos == "Grupo 1", "#FF0000",
ifelse(dato_mapa$Grupos == "Grupo 2", "#CCFF00",
ifelse(dato_mapa$Grupos == "Grupo 3", "#00FF66",
ifelse(dato_mapa$Grupos == "Grupo 4", "#0066FF",
ifelse(dato_mapa$Grupos == "Grupo 5","#CC00FF",0)))))
# Crear mapa
leaflet() %>% addPolygons(data = dato_mapa, opacity = 0.4, color = "#545454",weight = 1, fillColor = dato_mapa$colores,
fillOpacity = 0.4, label = ~NOMBRE_BAR,
highlightOptions = highlightOptions(color = "#262626", weight = 3, bringToFront = T, opacity = 1),
popup = paste("Barrio: ", dato_mapa$NOMBRE_BAR, "<br>", "Grupo: ", dato_mapa$Grupos)) %>%
addProviderTiles(providers$OpenStreetMap) %>%
addLegend(position = "bottomright", colors = thalia, labels = c("Grupo 1: Muy alta accidentalidad y mortalidad", "Grupo 2: Poca accidentalidad", "Grupo 3: Moderada accidentalidad", "Grupo 4: Alta accidentalidad", "Grupo 5: Muy poca accidentalidad"))Para Predecir la accidentalidad de la ciudad de medellín se ha dispuesto de la informacion suministrada por el IDEAM (Instituto de Hidrología, Meteorología y Estudios Ambientales), con esta se recaudo la medida de la pluviometría diaria y la humedad relativa a nivel del suelo en la ciudad de medellin.
A partir de esta información y la obtenida por el canal de datos abiertos de movilidad que publica la Alcaldía de Medellín en el portal GeoMedellín donde se exhiben los diferentes accidentes diarios y medidas del mismo,junto con la compilacion de las fechas festivas colombianas y generando un atributo para determinar si el accidente ocurrio en medio de horas pico, teniendo en cuenta que se desea predecir la accidentalidad por tipo de accidente a nivel semanal, mensual y diario, se procedio a separar esta base de datos en dos, de la historia de accidentalidad que comprende desde el 2014 hasta el 2017, y la del 2018, esta ultima para realizar validaciones de los modelos planteados, esta base de datos la llamamos accidentalidad_medellin_final.csv y puede descargarla de github.
Asi, se procedio a agrupar todos los accidentes sucedidos en un mismo Dia, o en una misma semana, o en un mismo mes,junto con el tipo de accidente (Choque, Atropello, Caida Ocupante, Volcamiento, Incendio, Otro) usando funciones de agregacion para los diferentes atributos expuestos anteriormente, puede ver el script y las bases de datos en aqui.
Ya agrupadas y listas las bases de datos por mes,semana y dia al igual que por tipo de accidentes, procederemos a explicar brevemente los modelos que se utilizaron para generar las predicciones del año 2018 y sus desempeños.
A continuación se dará una breve descripción de la manera en que se trataron y usaron las bases de datos para entrenar modelos predictivo sobre la cantidad de accidentes en un periodo de tiempo, para esto se evaluaron 5 tipos de modelos:
También es importante tener en cuenta las resoluciones temporales que se trabajarán los modelos y son las siguientes:
Los registros faltantes en las variables del clima se les realizó imputación basado en su media para la Humedad relativa y valor de cero para la Pluviometria. A continuación se presenta como queda conformado el conjunto de entrenamiento para la resolución semanal, de esta misma forma queda construido el nivel diario y mensual para los conjuntos de entrenamiento y validación.
require(lubridate)
require(tidyverse)
require(fastDummies)
require(xgboost)
require(Matrix)
require(gamlss)
require(skimr)
require(ranger)
library(FNN)
require(neuralnet)
require(gotop)
use_gotop()#==============================
# Imputación datos faltantes
#===============================
medellin$HUMEDAD_REL[is.na(medellin$HUMEDAD_REL)] <- mean(medellin$HUMEDAD_REL, na.rm = T)
medellin$PLUVIO[is.na(medellin$PLUVIO)] <- 0
#=======================================================
# Elección conjunto de entrenamiento y validación
#=======================================================
Train <- medellin[medellin$PERIODO %in% c(2014,2015,2016,2017),]
Test <- medellin[medellin$PERIODO == 2018,]
#=========================
# Puesta a Punto
#==========================
# Convertir ESPECIAL y HORA PICO en binarias
Train$ESPECIAL <- ifelse(Train$ESPECIAL=='SI',1,0)
Train$HORA_PICO <- ifelse(Train$HORA_PICO=='SI',1,0)
Test$ESPECIAL <- ifelse(Test$ESPECIAL=='SI',1,0)
Test$HORA_PICO <- ifelse(Test$HORA_PICO=='SI',1,0)
Train$CLASE <- as.factor(Train$CLASE)
levels(Train$CLASE) <- c('Atropello','Caida_ocupante', 'Choque','Choque',
'Incendio','Otro','Volcamiento')
Test$CLASE <- as.factor(Test$CLASE)
levels(Test$CLASE) <- c('Atropello','Caida_ocupante', 'Choque','Choque',
'Incendio','Otro','Volcamiento')
#====================================================
# Agrupamiento variables temporales para los modelos
#====================================================
# Grupo1: Mensual
Mensual <- Train %>% group_by(Mes = month(FECHA),CLASE) %>%
summarise(cantidad = length(FECHA),
festivos = sum(ESPECIAL),
Hora_Pico = sum(HORA_PICO),
ValorHumR = mean(HUMEDAD_REL, na.rm = T),
Pluvial = max(PLUVIO)) %>%
dummy_cols(select_columns = c('CLASE'),
remove_selected_columns = T)
MensualTs <- Test %>% group_by(Mes = month(FECHA),CLASE) %>%
summarise(cantidad = length(FECHA),
festivos = sum(ESPECIAL),
Hora_Pico = sum(HORA_PICO),
ValorHumR = mean(HUMEDAD_REL, na.rm = T),
Pluvial = max(PLUVIO)) %>%
dummy_cols(select_columns = c('CLASE'),
remove_selected_columns = T)
# Grupo 2: Diario
Diario <- Train %>% group_by(Dia = day(FECHA),CLASE) %>%
summarise(cantidad = length(FECHA),
festivos = sum(ESPECIAL),
Hora_Pico = sum(HORA_PICO),
ValorHumR = mean(HUMEDAD_REL, na.rm = T),
Pluvial = max(PLUVIO)) %>%
dummy_cols(select_columns = c('CLASE'),
remove_selected_columns = T)
DiarioTs <- Test %>% group_by(Dia = day(FECHA),CLASE) %>%
summarise(cantidad = length(FECHA),
festivos = sum(ESPECIAL),
Hora_Pico = sum(HORA_PICO),
ValorHumR = mean(HUMEDAD_REL, na.rm = T),
Pluvial = max(PLUVIO)) %>%
dummy_cols(select_columns = c('CLASE'),
remove_selected_columns = T)
# Grupo 3: Semanal
Semanal <- Train %>% group_by(Semana = week(FECHA),CLASE) %>%
summarise(cantidad = length(FECHA),
festivos = sum(ESPECIAL),
Hora_Pico = sum(HORA_PICO),
ValorHumR = mean(HUMEDAD_REL, na.rm = T),
Pluvial = max(PLUVIO)) %>%
dummy_cols(select_columns = c('CLASE'),
remove_selected_columns = T)
SemanalTs <- Test %>% group_by(Semana = week(FECHA),CLASE) %>%
summarise(cantidad = length(FECHA),
festivos = sum(ESPECIAL),
Hora_Pico = sum(HORA_PICO),
ValorHumR = mean(HUMEDAD_REL, na.rm = T),
Pluvial = max(PLUVIO)) %>%
dummy_cols(select_columns = c('CLASE'),
remove_selected_columns = T)
head(Semanal)Ahora, se describe el proceso uno a uno de los modelos utilizados y su implemantación en los datos:
XGBoost es una biblioteca optimizada de aumento de gradiente distribuida diseñada para ser altamente eficiente, flexible y portátil. Implementa algoritmos de aprendizaje automático en el marco de Gradient Boosting.
Boosting es una metodologia para mejorar la prediccion proporcionada por un arbol de decision. Como Bagging, Boosting es un procedimiento general que puede ser aplicado a muchos metodos de aprendizaje estadistico para regresion y clasicacion. En Boosting en vez de ajustar un arbol de decision grande, lo cual tiene el potencial de de sobre ajustar los datos, en forma secuencial se adicionan arboles donde cada uno es ajustado a los residuos del anterior.
Tome \(\hat{f}(x)=0\) y \(r_i=y_i\) para todo i en el conjunto de entrenamiento.
ahora, para un b=1,2,3,…,B se repiten los siguientes pasos:
Ajuste un arbol con d bifurcaciones (d+1 nodos terminales) a los datos (\(X,r\)).
Actualice \(\hat{f}\) adicionando una version ‘contraida’ del arbol nuevo, es decir:
\[\hat{f} \rightarrow \hat{f}+\lambda\hat{f^b}(x)\] - Actualice los residuos:
\[r_i=r_i-\lambda\hat{f^b}(x_i)\] Finalmente se obtiene un modelo Boosted
\[\hat{f}=\sum_{b=1}^B\lambda\hat{f^b}(x)\]Ahora, se genero un modelo XGboost que reducia la funcion de costo de la raiz del MSE, con un parametro de aprendizaje de \(\lambda=1\) y 100 arboles.
La forma en la que se entreno el modelo en R es como sigue:
#============================
# XGBoost
#============================
param1 <- list(booster = "gblinear", nthread = 100, eta = 0.1,
gamma = 0, max = 10, min_child_weight = 1,
max_delta_step = 2, lambda = 1,
alpha = 0, three_method = "exact",
objective = "reg:squarederror",
eval_metric = "rmse") # Parametros del Boosting
#-----------------
# XGboost: Mensual
#-----------------
dtrain <- xgb.DMatrix(as.matrix(Mensual[,-c(2)]),label = Mensual$cantidad)
dtest <- xgb.DMatrix(as.matrix(MensualTs[,-c(2)]),label = MensualTs$cantidad)
watchlist <- list(train=dtrain,test=dtest)
ModMesXGboost <- xgb.train(data = dtrain,nrounds = 200,
watchlist = watchlist,params = param1,
max.depth=2)
#------------------
# XGboost: Semanal
#------------------
dtrain <- xgb.DMatrix(as.matrix(Semanal[,-c(2)]),label = Semanal$cantidad)
dtest <- xgb.DMatrix(as.matrix(SemanalTs[,-c(2)]),label = SemanalTs$cantidad)
watchlist <- list(train=dtrain,test=dtest)
ModSemXGboost <- xgb.train(data = dtrain,nrounds = 100,
watchlist = watchlist,params = param1,
max.depth=2)
#-----------------
# XGboost: Diario
#-----------------
dtrain <- xgb.DMatrix(as.matrix(Diario[,-c(2)]),label = Diario$cantidad)
dtest <- xgb.DMatrix(as.matrix(DiarioTs[,-c(2)]),label = DiarioTs$cantidad)
watchlist <- list(train=dtrain,test=dtest)
ModDiaXGboost <- xgb.train(data = dtrain,nrounds = 100,
watchlist = watchlist,params = param1,
max.depth=2)
#------------------------
# Validacion MSE XGboost
#------------------------
require(Matrix)
MSEXGboostMes <- mean((predict(ModMesXGboost,
as(as.matrix(MensualTs[,-2]),"dgCMatrix"))-MensualTs$cantidad)^2)
MSEXGboostSem <- mean((predict(ModSemXGboost,
as(as.matrix(SemanalTs[,-2]),"dgCMatrix"))-SemanalTs$cantidad)^2)
MSEXGboostDia <- mean((predict(ModDiaXGboost,
as(as.matrix(DiarioTs[,-2]),"dgCMatrix"))-DiarioTs$cantidad)^2)finalmente se entreno el modelo para las bases de agrupadas por Mes, Semana y Dia asi como por tipo de accidente, se obtuvieron los siguientes resultados en la validacion cruzada por el metodo de conjunto de validacion sobre el año 2018.
Los GAMLSS (Generalized Additive Models for Location,Scale and Shape) son modelos de regresión distributiva univariante, donde todos los parámetros de la distribución asumida para la respuesta se pueden modelar como funciones aditivas de las variables explicativas.
Es una mejora del modelo GAM (Generalized Additive Model), definido asi:
\(Y_i = \beta_0+f_1(X_{i1})+f_2(X_{i2})+...+f_p(X_{ip})+\epsilon\)
Para el modelo GAMLSS consideramos los siguientes vectores de parametros:
Asumiendo que Y es la variable respuesta y \(D\) es alguna distribucion, el modelo GAMlSS se plantea como:
\[Y \sim D(\mu,\sigma,\nu,\tau)\] \[g_1(\mu) = \beta_{01}+f_{11}(X_1)+f_{12}(X_2)+ ... + f_{1p}(X_p)\] \[g_2(\mu) = \beta_{02}+f_{21}(X_1)+f_{22}(X_2)+ ... + f_{2p}(X_p)\] \[g_3(\mu) = \beta_{03}+f_{31}(X_1)+f_{32}(X_2)+ ... + f_{3p}(X_p)\] \[g_4(\mu) = \beta_{04}+f_{41}(X_1)+f_{42}(X_2)+ ... + f_{4p}(X_p)\]Con lo que a diferentes funciones \(f\) o tambien llamadas links y una distribucion \(D\) podremos generar un modelo que prediga la cantidad de accidentes en un tiempo determinado (experimento \(Poisson\))
Para nuestro caso particular utilizamos como funcion \(D\) una \(Poisson\) y el paquete gamlss posee una funcion de ajuste aproximado para \(D\) fitDist() con la que tambien se modelo una \(Exponential\), siendo esta ultima la que mejor desempeño en prueba obtuvo con los siguientes valores por el metodo de validacion en el año 2018.
La forma en la que se entreno el modelo en R es como sigue:
#=========================================================================
# Regresión lineal generalizada usando GAMLSS
#==========================================================================
#------------------
# RLGamlss: Mensual
#-------------------
ModMenGam <- gamlss(formula = cantidad~.,
sigma.formula = ~.,
data = Mensual,
family = EXP)
#-------------------
# RLGamlss: Semanal
#--------------------
ModSemGam <- gamlss(formula = cantidad~.,
sigma.formula = ~.,
data = na.omit(Semanal),
family = EXP)
#--------------------
# RLGamlss: Diario
#--------------------
ModDiaGam <- gamlss(formula = cantidad~.,
sigma.formula = ~.,
data = na.omit(Diario),
family = EXP)
#-------------------------------
# Validacion MSE para RLGamlss
#-------------------------------
MSEgamlssMes <- mean((predictAll(ModMenGam, MensualTs)$mu-MensualTs$cantidad)^2)
MSEgamlssSem <- mean((predictAll(ModSemGam, SemanalTs)$mu-SemanalTs$cantidad)^2)
MSEgamlssDia <- mean((predictAll(ModDiaGam, DiarioTs)$mu-DiarioTs$cantidad)^2)Es un metodo de ensamble de arboles de clasificacion o regresion, donde se consideran arboles sin poda cada arbol es entrenado sobre una muestra boostrap de la base de interes, y a su vez en cada split del arbol se utiliza una muestra aleatoria de los atributos en la base de interes.
En nuestro modelamiento se utilizaron 350 arboles de regresion y muestras de 4 variables en cada split de la cantidad de accidentes en cierto periodo respecto a las demas covariables, ademas se opto por una libreria considerablemente mas optimizada ranger lo que permite realizar muchas pruebas de bosques con diferentes valores para los parametros en poco tiempo.
asi, nuestro bosque aleatorio, fue entrenado en base a datos del 2014 al 2017 y validado con datos del año 2018, obteniendo las metricas del MSE por conjunto de validacion del año 2018 como sigue.
La forma en la que se entreno el modelo en R es como sigue:
#====================================
# Random Forest
#====================================
#------------------------
# Random Forest: Mensual
#------------------------
ModMenRanger <- ranger(cantidad~.,
data = na.omit(Mensual),
num.trees = 350, mtry = 4)
predRFMes <- predict(ModMenRanger, na.omit(MensualTs))
#-----------------------
# Random Forest: Semanal
#-----------------------
ModMenRanger <- ranger(cantidad~.,
data = na.omit(Semanal),
num.trees = 350,mtry = 4)
predRFSem <- predict(ModMenRanger, na.omit(SemanalTs))
#-----------------------
# Random Forest: Diario
#-----------------------
ModMenRanger <- ranger(cantidad~.,
data = na.omit(Diario),
num.trees = 350,mtry = 4)
predRFDia <- predict(ModMenRanger, DiarioTs)
#---------------------------------
# Validacion MSE Random Forest
#---------------------------------
MSERangerMes <- mean(((predRFMes$predictions)-MensualTs$cantidad)^2)
MSERangerSem <- mean(((predRFSem$predictions)-SemanalTs$cantidad)^2)
MSERangerDia <- mean(((predRFDia$predictions)-DiarioTs$cantidad)^2)Este metodo permite obtener el valor de una prediccion a travez de la media de sus k vecinos mas cercanos (principalmente metricas euclidianas), es decir, para hallar el valor predicho de una observacion \(X_i\) en \(y_i\) y dado un valor de k, se promedian los valores \(y\) de las k observaciones mas cercanas a la de interes, este promedio es el valor predicho para esta observacion.
En nuestro caso se probaron valores de k = 1,2,3,…,50 para cada base de datos agrupada, asi se encontro el K-optimo y con este se obtuvo los siguientes resultados de validacion por el metodo de conjunto de validacion con datos del 2018.
La forma en la que se entreno el modelo en R es como sigue:
#========================================
# K-Nearest Neighbor Regression
#========================================
#-----------------
# KNN: Mensual
#-----------------
MSEKnnMes <- vector()
for(i in 1:50){
ModMesKnn <- knn.reg(train = na.omit(Mensual[,-c(2,3)]),
test = na.omit(MensualTs[,-c(2,3)]),
y = Mensual$cantidad,k=i)
MSEKnnMes[i] <- mean(((ModMesKnn$pred)-na.omit(MensualTs)$cantidad)^2)
}
plot(1:50,MSEKnnMes, type = 'l')
paste("el valor de K que minimisa el MSE es:",which.min(MSEKnnMes))
#-------------
# KNN: Semanal
#-------------
MSEKnnSem <- vector()
for(i in 1:50){
ModSemKnn <- knn.reg(train = na.omit(Semanal[,-c(2,3)]),
test = na.omit(SemanalTs[,-c(2,3)]),
y = Semanal$cantidad,k=i)
MSEKnnSem[i] <- mean(((ModSemKnn$pred)-na.omit(SemanalTs)$cantidad)^2)
}
plot(1:50,MSEKnnSem, type = 'l')
paste("el valor de K que minimisa el MSE es:",which.min(MSEKnnSem))
#-------------
# KNN: Diario
#-------------
MSEKnnDia <- vector()
for(i in 1:50){
ModDiaKnn <- knn.reg(train = na.omit(Diario[,-c(2,3)]),
test = na.omit(DiarioTs[,-c(2,3)]),
y = Diario$cantidad,k=i)
MSEKnnDia[i] <- mean(((ModDiaKnn$pred)-na.omit(DiarioTs)$cantidad)^2)
}
plot(1:50,MSEKnnDia, type = 'l')
paste("el valor de K que minimisa el MSE es:",which.min(MSEKnnDia))Las redes neuronales son un modelo computacional vagamente inspirado en el comportamiento observado en su homólogo biológico. Consiste en un conjunto de unidades, llamadas neuronas artificiales, conectadas entre sí para transmitirse señales. La información de entrada atraviesa la red neuronal (donde se somete a diversas operaciones) produciendo unos valores de salida.
El algoritmo basicamente genera en principio unos pesos \(W_1\) aleatorios entre las capas definidas por el usuario. De forma resumida, estos pesos se multiplican por cada registro y este resultado se le aplica una funcion no lineal a la que llamaremos activacion, y con esto obtendremos los valores de la primer capa de la red, a estos valores se le asignan unos pesos \(W_2\) aleatorios para realizar de nuevo la operacion anteriormente descrita y obtener la segunda capa de la red asi sucesivamente entre todas las capas de la estructura definida por el usuario.
Cuando se llega a la ultima capa llamada capa de salida(definida por el usuario) y despues de haber calculado todas las neuronas y sus activaciones (llamado propagacion hacia adelante), el algoritmo calcula una funcion de costo o funcion objetivo, con la intencion de reducir esta funcion (generalmente una distancia ponderada de el valor de la capa de salida a un vector \(y\) llamado Target) se recurre a recalcular los \(W\) para que se minimice la funcion de costo, a esto se le llama propagacion hacia atras, cabe aclarar que en cada activacion de las neuronas se generan vectores \(b\) que se suman a cada neurona.
En nuestro caso particular y notando que el metodo de redes es computacionalmente demandante, generamos las redes neuronales que se podian correr en nuestra maquina, es decir utilizamos una arquitectura de red tal que pudiese correrse en nuestra maquina y fueron las siguientes:
Con estas calculamos los MSE de validacion para el periodo 2018, obteniendo:
A continuación verá un resumen gráfico de las validaciones de los modelos teniendo como medida la Raiz del MSE. Se resalta la notoriedad de los modelos KNN y XGboost en la reducción del MSE sobre todos los demás usando como validación el 2018, siendo el XGboost el que logra mejor desempeño para los 3 periodos, por lo tanto este será el modelo para realizar las predicciones.
Modelo <- c("XGBoost","XGBoost","XGBoost",
"GAMLSS","GAMLSS","GAMLSS",
"RandomForest","RandomForest","RandomForest",
"KNN","KNN","KNN",
"RedNeuronal","RedNeuronal","RedNeuronal")
Periodo <- c("Mensual","Semanal","Diario",
"Mensual","Semanal","Diario",
"Mensual","Semanal","Diario",
"Mensual","Semanal","Diario",
"Mensual","Semanal","Diario")
MSE <- c(7882.709,507.5032,1522.026,
1151768,22477.15,96030.66,
1404420,35593.04,301954.1,
98177.22,4888.43,16135.33,
11520464,0,1356718)
Evaluacion <- data.frame(Modelo = Modelo, Periodo = Periodo, MSE = MSE)
ggplot(data = Evaluacion, aes(x = Modelo, y = sqrt(MSE), fill = Periodo)) +
geom_bar(stat = "identity", position = "dodge") +
theme_minimal() +
xlab("Modelo") +
ylab("RMSE")Dado que en los anteriores modelos se obtiene con los accidentes agrupados por dias, semana y mes, se necesita un modelo que relacione una a uno los dias teniendo en cuenta las fechas. Dado que los accidentes representa un conteo, se entrenará un modelo de regresión generalizada usando la distribución Poisson, para las variables FECHA, CLASE, DIA_NOMBRE, PERIODO, ESPECIAL.
#=============================
# MODELO GENERALIZADO POISSON
#=============================
Train$CLASE <- toupper(Train$CLASE)
Test$CLASE <- toupper(Test$CLASE)
grupo_1 <- Train %>% group_by(FECHA, CLASE, DIA_NOMBRE, PERIODO, ESPECIAL) %>%
summarise(Accidentes = n()) %>%
arrange(FECHA)
# Modelo Poisson
modelo_poisson <- glm(Accidentes~FECHA+CLASE+DIA_NOMBRE+PERIODO+ESPECIAL,
data = grupo_1, family = "poisson")
summary(modelo_poisson)##
## Call:
## glm(formula = Accidentes ~ FECHA + CLASE + DIA_NOMBRE + PERIODO +
## ESPECIAL, family = "poisson", data = grupo_1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -8.0990 -1.0024 -0.0941 0.7931 6.0641
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.199e+02 1.722e+01 6.960 3.39e-12 ***
## FECHA 1.877e-04 2.315e-05 8.110 5.06e-16 ***
## CLASECAIDA_OCUPANTE -1.345e-01 1.122e-02 -11.987 < 2e-16 ***
## CLASECHOQUE 1.890e+00 8.221e-03 229.948 < 2e-16 ***
## CLASEINCENDIO -2.465e+00 2.427e-01 -10.159 < 2e-16 ***
## CLASEOTRO 6.831e-02 1.066e-02 6.411 1.45e-10 ***
## CLASEVOLCAMIENTO -1.093e+00 1.558e-02 -70.175 < 2e-16 ***
## DIA_NOMBREJUEVES 4.584e-01 1.007e-02 45.501 < 2e-16 ***
## DIA_NOMBRELUNES 4.244e-01 1.022e-02 41.505 < 2e-16 ***
## DIA_NOMBREMARTES 4.882e-01 1.001e-02 48.779 < 2e-16 ***
## DIA_NOMBREMIERCOLES 4.636e-01 1.005e-02 46.134 < 2e-16 ***
## DIA_NOMBRESABADO 3.979e-01 1.017e-02 39.125 < 2e-16 ***
## DIA_NOMBREVIERNES 5.057e-01 9.973e-03 50.700 < 2e-16 ***
## PERIODO -6.002e-02 8.733e-03 -6.873 6.28e-12 ***
## ESPECIAL -2.453e-01 1.020e-02 -24.060 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 199129 on 7246 degrees of freedom
## Residual deviance: 13955 on 7232 degrees of freedom
## AIC: 45706
##
## Number of Fisher Scoring iterations: 4
Para la evaluacion del modelo se usa nuevamente el MSE y evaluamos para el conjunto de entrenamiento y validación y calculamos la variación porcentual entre los dos:
#MSE para entrenamiento
prediccion_Train <- predict(modelo_poisson, type = "response")
MSE_Train <- mean((prediccion_Train - grupo_1$Accidentes)^2)
# MSE Validacion
grupo_validacion <- Test %>% group_by(FECHA, CLASE, DIA_NOMBRE,
PERIODO, ESPECIAL) %>%
summarise(Accidentes = n()) %>%
arrange(FECHA)
prediccion_Test <- predict(object = modelo_poisson, newdata = grupo_validacion[,-6],
type = "response")
MSE_Test <- mean((prediccion_Test - grupo_validacion$Accidentes)^2)
print(c("MSE Test" = MSE_Test, "MSE Train" = MSE_Train))## MSE Test MSE Train
## 55.69298 57.50773
## Variacion porcentual
## 1.629249
Dado que los modelos son temporales por lo que entre más alejado sea el horizonte de pronóstico menos acertada será la predicción se tomará el año 2019 para predecir con los modelos entrenados con 4 años (2014-2017). Los datos pronosticados se guardarán en un csv que será utilizado para la visualización e interacción del usuario en la aplicación VEHAPP
#====================================
# Prediccion 2019 con modelo poisson
#=====================================
Datos_2019 <- read.csv("Datos/Otros/Predecir_2019.csv", sep = ",", encoding = "UTF-8")
Datos_2019 <- Datos_2019[,-1]
Datos_2019$CLASE <- gsub(pattern = "CAIDA OCUPANTE", replacement = "CAIDA_OCUPANTE", x = Datos_2019$CLASE)
Datos_2019$FECHA <- as.Date(Datos_2019$FECHA, format = "%Y-%m-%d")
Datos_2019$CLASE <- as.factor(Datos_2019$CLASE)
Datos_2019$DIA_NOMBRE <- as.factor(Datos_2019$DIA_NOMBRE)
Datos_2019$PERIODO <- as.integer(Datos_2019$PERIODO)
prediccion2019 <- predict(object = modelo_poisson, newdata = Datos_2019,
type = "response")
PREDICCIONES_DIARIAS_2019 <- Datos_2019 %>%
mutate(ACCIDENTES = round(prediccion2019,0))
#write.csv(PREDICCIONES_DIARIAS_2019, "Predicciones_diarias_2019.csv")En la caracterización de los accidentes por barrios se encontró que la diferencia más significativa entre estos es la cantidad de accidentes que ocurren y se llegó a la conclusión que entre más cercano se está a la zona centro más propenso está a sufrir un accidente y que las vias rapidas de la ciudad son centro constante de accidentes también.
Como observación final, el uso de técnicas estadísticas depende mucho de la calidad del dato y de como estos se hayan recolectado y tratado.
– Software R